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Abstract 

Using a multi-scaled, chaotic flow known as the KS model of turbulence [7] , we investigate the dependence of Lyapunov 
exponents on various characteristics of the flow. We show that the KS model yields a power law relation between 
the Reynolds number and the maximum Lyapunov exponent, which is similar to that for a turbulent flow with the 
same energy spectrum. Our results show that the Lyapunov exponents are sensitive to the advection of small eddies 
by large eddies, which can be explained by considering the Lagrangian correlation time of the smallest scales. We 
also relate the number of stagnation points within a flow to the maximum Lyapunov exponent, and suggest a linear 
dependence between the two characteristics. 
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1. Introduction 

Measures of stretching, such as the Lyapunov ex- 
ponent, are important tools for understanding the 
nature of dynamical systems. For example, the max- 
imum Lyapunov exponent can provide information 
about the complexity of an attractor via the Kaplan- 
Yorke dimension [6] , or the rate of loss of information 
in the system [16]. The use of Lyapunov exponents 
in turbulent flows is far too great to list here; exam- 
ples range from probing the onset of turbulence [3] , 
to detecting inhomogeneity in hydromagnetic con- 
vection [9] . Another interesting application arises in 
dynamo theory; it was shown [15] that Lyapunov ex- 
ponents provide an upper bound for the growth rate 
of a fast dynamo, and also a non-trivial combination 
of Lyapunov exponents gives an exact growth rate 
for the small scale turbulent dynamo [5] . 
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In this work we use a model turbulent flow, known 
as the Kinematic Simulation (KS) model, that has 
been primarily used as a Lagrangian model of tur- 
bulence [7,8]. An important feature of KS is that 
it allows full control of the energy spectrum; more- 
over, its simple analytic structure means that nu- 
merical differentiation is not required in calculating 
the Lyapunov exponents. The KS model has been 
shown to be in good agreement with results obtained 
from direct numerical simulations (DNS) of turbu- 
lent flows, particularly with respect to Lagrangian 
statistics such as two-particle dispersion [8,10,11]. 
The use of the model is spreading rapidly to many 
other areas such as aeroacoustics and biomechanics. 
This flow has also been shown to be a hydromagnetic 
dynamo [17]. Motivated by the success of the KS 
model and its applications in magnetohydrodynam- 
ics, our aim is to check the agreement between the 
model and turbulent flows with respect to Lyapunov 
exponents. It has been shown that the KS model 
exhibits Lagrangian chaos [8,10]; we shall quantify 
this feature using the largest Lyapunov exponent. 
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2. The velocity field 

The KS model prescribes the flow velocity at a 
position X and time t through the summation of 
Fourier modes with randomly chosen parameters. 
These modes are mutually independent, therefore 
the advection of small eddies by large eddies is not 
included in the model. More precisely, the velocity 
field is prescribed to be [11] 

N 



Since the maximum value of E{k) and the inte- 
gral scale remain unchanged in the models discussed 
here, any variation in Re can be thought to be caused 
by changes in the fluid viscosity. Fig. 1 shows the 
energy spectrum of the KS flow, obtained numer- 
ically after fast Fourier transforming u calculated 
from Eq. (1) on a 128^ mesh. We also show a slice, 
in the z plane, of the corresponding vorticity field, 
with velocity vectors. 



u(x, i) = ^ ( A„ X k„ cos ipn+BnX k„ sin V)™) , 3. The Lyapunov exponents 



n=l 

(1) 

where ijjn = k„ • x + ujnt and N is the number of 
modes. The unit vectors k„ are chosen randomly, 
and k„ = fc„k„ where fc„ is the wavenumber of the 
^th ]-Qode. The construction of A„ and B„, which 
are time independent, is explained in the appendix. 
Even though the parameters of the flow are chosen 
randomly, they do not necessarily change with time, 
so the flow is not necessarily random. We adopt a 
normalised energy spectrum of the KS flow E{k), 
which is a modification of the von Karman energy 
spectrum, 

E{k) = k*{l + fc2)-(2+P/2)e-i/2(fc/fc«)^^ (2) 

which reduces to E{k) oc k~^ in the inertial range 
1 ^ fc <C /cat, with /c = 1 at the integral scale; p = 
5/3 produces the Kolmogorov spectrum. As men- 
tioned previously, a useful feature of the KS model 
is the ability to vary the slope p in the inertial range. 
The flow is incompressible and time dependent; the 
frequency of the n**^ mode, a;„ is inversely propor- 
tional to its turnover time, 

a;„ = ^klE{kn). (3) 

It is convenient to write the unit vector k„ as 



\/l - CnCos6'„ 
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where, 6'„ G [0, 27r) and C„ G [—1,1], are uniformly 
distributed random numbers, to ensure that k„ are 
isotropically distributed. With 
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the effective Reynolds number is introduced us- 
ing the requirement that the dissipation and eddy 
turnover times are equal to each other at A; = fcjv, 

Re = {kN/k^)^P+'^y^. (6) 



To obtain the spectrum of Lyapunov exponents, 
Aj, we measure the average rates of exponential di- 
vergence of nearby fluid particle trajectories. If the 
system is chaotic, at least one Lyapunov exponent is 
positive. The procedure to calculate the Lyapunov 
exponents consists of monitoring the evolution of an 
infinitesimal fluid sphere moving with the flow. The 
sphere, deformed by the flow, rapidly becomes an 
ellipsoid. Then the Lyapunov exponents are defined 
as 

1 - P^{t) 
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where Pi{t) is the ellipsoid's i'^ principal axis, and 
i = 1,2,3. Another feature of the KS flow is that 
it is time reversible (unlike 'real' turbulence), there- 
fore the second Lyapunov exponent vanishes [2,1]. 
We now consider two remaining exponents, which 
must have opposite signs, Ai = —A3, since the flow 
is incompressible (V • u = 0), the sum of the Lya- 
punov exponents must be zero. Hence we only need 
to calculate one exponent, A = max(Ai). Following 
Wolf et al. [18], consider a sphere whose centre, at 
xq, moves along a trajectory defined by 



— Xo(t) = u(xo,t). 



(8) 



with u obtained from Eq. (1). As the sphere follows 
a trajectory in the fiow, its shape is deformed to an 
ellipsoid by stretching and compression. To the lin- 
ear approximation in the sphere radius rj, xq remains 
the centre of the deformed ellipsoid. Positions of the 
points on the surface of the sphere r] = x — Xo, where 
Xo is the position of the centre of the ellipsoid, obey 
the linearised equations of motion 



dt 



rii{t) = Dijr]j, 



(9) 



where 



dui/dxj and the summation conven- 
tion is assumed. We integrate Eq. (8) and (9) nu- 
merically, normalising rj at regular intervals as to 
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Fig. 1. (a) The energy spectrum, E{k), showing the imposed p = 5/3 slope, as obtained by Fourier transform of Eq. (1) with 
N = 20, fei = 10 and fcjv = 400. (b) Slice in the z plane of the vorticity field generated by taking curl of the velocity field 
from (a), lighter shading indicates higher vorticity. Velocity vectors axe shown in white. 
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Fig. 2. The average Lyapunov exponent obtained from av- 
eraging over 500 particles within several realizations of the 
flow (1) is shown with dashed lines. The average over 50 re- 
alisations of the flow is shown with solid curve. The flow has 
AT = 20, fei = 6 and fejv = 90. 



keep the linearisation valid. We then take the tem- 
poral average of the magnitude of rj to recover the 
maximum Lyapunov exponent. Finally, wc average 
the results over 500 particles to improve statistics. 
Since detailed behavior of A can vary significantly 
between different realisations of the flow, wc further 
take an ensemble average over 50 diflFerent realisa- 
tions of the KS model with the same non-random 
parameters. The results of one such run are shown 
in Fig. 2. Before beginning the simulations, the code 
was tested by computing Lyapunov exponents for 
some well known chaotic flows [13] . 



4. Scaling of the Lyapunov exponent with 
the Reynolds number 

We now determine how A scales with the proper- 
ties of the flow. We begin by considering the relation- 
ship between the maximum Lyapunov exponent. A, 
and the Reynolds number, Re. Following Ruelle [12], 
we increase Re by introducing smaller scales (i.e., 
increasing fcjv), keeping the same number of modes 
and the same fci = 1. In this way we introduce mo- 
tions with higher velocity shear rate. The maximum 
Lyapunov exponent is the modulus of the long-term 
average of the velocity gradient in the Lagrangian 
frame. In fully developed turbulence this is related 
to the turnover time of the smallest eddies, where 
the turnover time of an eddy of size £ is 



T{i) 



1-h 



(10) 



where L is the integral scale, U the corresponding 
velocity (with tl = L/U) and h the Holder exponent 
of the velocity field, introduced as 
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Then A scales with Re as 
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where a = {1 — h)/{l + h) and irj = 27r/A:jv is the 
Kolmogorov length scale. For the Kolmogorov spec- 
trum, p = 5/3, we have h = 1/3 and a = 1/2. 
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Fig. 3. (a) A against Re with error bars, with A'^ = 40 in each reahzation and p = -5/3. Errors arc calculated from the scatter 
in A between different trajectories in each realisation, and between different realisations, (b) As in (a), but with the line of 
best fit shown dashed, giving a » 0.38. 



Our numerical simulations, illustrated in Fig. 3, 
show that the largest Lyapunov exponent of the KS 
model scales as A oc Rc" "^*, thus a is smaller than 
Ruelle's prediction. This difference may arise from 
the lack of sweeping of the small eddies by the large 
eddies in the KS model. In KS, like in real turbu- 
lence, velocity is determined by the large scale 'ed- 
dies', and velocity gradients are determined by the 
small scale motions. But, the Lyapunov exponent 
is related to velocity gradients in the Lagrangian 
frame, hence the lack of advection of the small eddies 
is important. Because of the lack of advection, the 
Lagrangian correlation time of eddies of size l//c„ is 
1/kn, not l/uJn ~ The maximum Lyapunov 

exponent is related to this correlation time, hence 
there is a reduction in the expected value of A, that 
grows with Re. We are grateful to an anonymous 
referee who suggested this explanation. 

Hence we must be careful when applying the KS 
model to area's where stretching is important, es- 
pecially if we are comparing results with DNS, and 
scaling with Re. Dynamo action is one of such areas. 

5. Lyapunov exponents, spectral slope and 
stagnation points 

Since the shear rate increases with the wave num- 
ber for a sufficiently steep spectrum, the main cause 
of stretching are the small scale motions, described 
well by the KS model. Therefore, we should not be 
too worried about the discrepancy described at the 
end of the previous section. Of our primary inter- 
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Fig. 4. The maximum Lyapunov exponent A plotted against 
the spectral slope p. 

est here (motivated by the saturation of the fluctu- 
ation dynamo) is the effect of the slope of the en- 
ergy spectrum plays on stretching. Cattaneo et al. [4] 
showed, using a simple chaotic velocity field, that as 
a dynamo saturates, Lagrangian chaos in the flow is 
suppresscid. With a complex multi-scaled flow such 
as turbulence this effect would remove energy from 
small scales first, stimulating us to study the effect 
of steepening the spectrum on A. To study how A 
depends on p, we fix k„ and kn/ki and change the 
spectral slope p in Eq. (2). The resulting values of 
A, shown in Fig. 4, have been obtained by averag- 
ing over 500 particles in the same fiow. These results 
confirm that the largest Lyapunov exponent is con- 
trolled mainly by the smaller scales: as less energy is 
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Fig. 5. Our estimate of the scaling parameter a plotted 
against the spectral slope p. The dashed line shows predicted 
values of a based on Ruelle's work. 

given to the small scales, A decreases. The next logi- 
cal step is to investigate the effect of p on the scaling 
parameter a. Indeed, h = ^{p — 1) in Eq. (11) and 
then 

AocRe(3-f)/(i+f\ (13) 

from Eq. (12). This dependence is shown dashed in 
Fig. 5 along with the results from our simulations. 
As before, our estimates of a are smaller, but we ex- 
pect this due to the absence of sweeping. If KS is to 
be used in areas sensitive to stretching, we suggest 
the inclusion of the advcc;tion of small eddies [7] . We 
would then expect the value of a to be closer to Ru- 
elle's prediction. It can be expected that the veloc- 
ity shear rate, and hence the amount of stretching, 
is maximum near stagnation points where velocity 
changes most rapidly in space. The number of stag- 
nation points per unit volume (and hence the total 
magnitude of velocity shear) is sensitive to the num- 
ber of modes in the KS model. The smaller is the 
value of kn+i — kn for large n, the more numerous 
will be the stagnation points [14], see Fig. 6. Al- 
though our velocity field is time-dependent, it has 
infinite correlation time, so that the density of stag- 
nation points varies little in time. We computed the 
maximum Lyapunov exponent for fixed Un /ki and 
p, but with an increasing number of modes, TV. 
Since the number of stagnation points in the flow 
increases with N, we expect that A increases too. 
Fig. 7 confirms this and also shows that A saturates 
as N grows. The abundance of stagnation points in 
the flow can be quantified using what we call the 
volume filling factor fv{(-r]) calculated as 



Fig. 6. Ux {x) shown for 100 modes in the solid line and and 10 
modes in the dashed line. With an increased number of modes 
the probability of each component being zero at a point in- 
creases, hence an increase in the number of stagnation points. 
A diagnostic of the shear rate, S = \ [ (summation 

convention assumed) also grows with A'': S = 11. 26 ±1.29 
and 17.75 ± 0.44, for N=10 and 100, respectively (averaged 
over 500 realisations). 
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Fig. 7. Plot of A against N , showing saturation of A as AT 
increases 

fv{(^v) = ]y3 XI -H"(0-05Mrms " \u{Xi, Vj , Zk)\) , 

(14) 

where the summation is extended over all the mesh 
points in the computational box, Nf, is the number 
of mesh points along each direction, Uims is the root- 
mean-square velocity, and H is the Heaviside func- 
tion, 

{0, ifa;<0, 
(15) 
1, ifa;>0. 
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Fig. 8. Error bar plot of fvi^r) against A with line of best 
fit, A Ri 63.5/y(£T) — 7.7, suggesting a linear relationship. 
Calculating the error in A is explained in Fig. 3, for error in 
fy we take error between realisations. 



Thus defined, this the fractional volume of the re- 
gion where |u| < O.OSurms- Wc calculated fvi^ri) in 
a KS flow with 10 < < 32 modes with iV^ = 128 
in a box whose size is in/ki, the largest scale in the 
flow, and wc average over 500 realisations of the flow. 
Fig. 8 shows a scatter plot of A versus for an 

increasing number of modes in the model. Despite 
the large error bars, the data suggests a linear rela- 
tionship. The magnitude of the errors may be due 
to the variation between different realisations of the 
flow. 



6. Conclusions 

Our simulations confirm that the KS model re- 
produces reasonably well the stretching properties 
of turbulent flows, including the scaling of the max- 
imum Lyapunov exponent with the Reynolds num- 
ber despite the fact that the latter can only be intro- 
duced formally in the KS model. The model values of 
the Lyapunov exponent are always smaller than the- 
oretical predictions for turbulent flows: we attribute 
this to the lack of advection of the small eddies by 
the large eddies in the model. Hence we must be 
careful when applying the KS model to area's where 
stretching is important. Finally our results suggest 
a linear relationship between the number of stagna- 
tion points in the flow and A. 



We choose A„, and B„ randomly, imposing or- 
thogonality with k„, which gives the required spec- 
trum as 

|A„xk„|=^„, (A.l) 

we proceed in the same fashion for B„. We then 
choose 



Bn 



2E{k„)Akn 



(A.2) 



This ensures 



tT / 7^\^\^dV= / E{k)dk^y^E{kn)Ak 

y Jv^ Jo 



where AA;„ is given by 

k2 - ki 



Akn 



2 ' 

kn - kn-i 



n = 1, 



l<n<N, 



n = N. 



(A.3) 



(A.4) 
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